rm(list = ls())
here::i_am(file.path("code", "10_war_experience.R"))
library(here)
source(here("code", "config.R"))

est_df <-  read_parquet(here("data", "analysis", "analysis.parquet"))

est_df <- est_df |>
  mutate(
    is_combat_div = (combat_div_ships == 92) | (combat_div_ships == 93),
    not_overseas = (vet_combined == TRUE) & (vet_ships == FALSE),
    overseas_noncombat = (vet_ships == TRUE) & (combat_div_ships == 0),
    div92 = (vet_ships == TRUE) & (combat_div_ships == 92),
    div93 = (vet_ships == TRUE) & (combat_div_ships == 93)
  )

p <- est_df |>
  mutate(
    grp = case_when(
      (vet_combined == TRUE) & is.na(combat_div_ships) ~ "not_overseas",
      (vet_combined == TRUE) & (combat_div_ships == 0) ~ "noncombat",
      (vet_combined == TRUE) & (combat_div_ships == 92) ~ "92",
      (vet_combined == TRUE) & (combat_div_ships == 93) ~ "93"
    )
  ) |>
  filter(!is.na(grp)) |>
  group_by(grp) |>
  summarize(
    est = mean(is_naacp),
    n = n(),
    se = sqrt(var(is_naacp)/n)
  ) |>
  mutate(
    grp = factor(
      grp,
      levels = c(
        "not_overseas",
        "noncombat",
        "92",
        "93"
      ),
      labels = c(
        "Noncombat,\nnot overseas",
        "Noncombat,\noverseas",
        "Combat,\n92nd Div.",
        "Combat,\n93rd Div."
      )
    )
  ) |>
  ggplot(aes(x = grp, y = est, ymin = est - 1.96 * se, ymax = est + 1.96 * se)) +
  my_theme() +
  labs(
    x = "",
    y = "Proportion in NAACP"
  )

p_color <- p + 
  geom_col(width = .5, fill = "#F8766D") +
  geom_errorbar(width = .1, color = "#333333")
ggsave(file.path(fig_dir, "color", "naacp_wartime_exp.pdf"), plot = p_color, width = 5, height = 3)

p_bw <- p + 
  geom_col(width = .5) +
  geom_errorbar(width = .1)
ggsave(file.path(fig_dir, "bw", "naacp_wartime_exp.pdf"), plot = p_bw, width = 5, height = 3)

vet_types <- c("not_overseas", "overseas_noncombat", "div92", "div93")

res <- feols(
  is_naacp ~ sw(not_overseas, overseas_noncombat, div92, div93) |
    sw0(birthyr_cards + exemption^married_cards + farmer + laborer + farm_laborer),
  data = est_df,
  vcov = "hetero"
)

p_color <- res |>
  map_dfr(broom::tidy, .id = "spec") |>
  filter(term != "(Intercept)") |>
  left_join(
    sapply(vet_types, \(x) est_df |> filter(get(x) == FALSE) |> 
             summarize(mean(is_naacp)) |>
             as.numeric() ) |>
      as_tibble(rownames = "term") |>
      mutate(term = paste0(term, "TRUE")),
    by = "term"
  ) |>
  mutate(
    controls = factor(
      ifelse(grepl("fixef: 1", spec), "No", "Yes"),
      levels = c("No", "Yes")
    ),
    estimate = value + estimate,
    conf.high = estimate + 1.96 * std.error,
    conf.low = estimate - 1.96 * std.error,
    term = factor(
      term,
      levels = c("not_overseasTRUE", "overseas_noncombatTRUE", "div92TRUE","div93TRUE"),
      labels = c("Noncombat,\nnot overseas",
                 "Noncombat,\noverseas",
                 "Combat,\n92nd Div.",
                 "Combat,\n93rd Div."))
  ) |>
  ggplot(aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high, fill = controls)) +
  geom_col(position = position_dodge(.75), width = .7) +
  geom_errorbar(position = position_dodge(.75), width = .1, color = "black") +
  my_theme() +
  labs(
    x = "",
    y = "Proportion in NAACP",
    fill = "Controls"
  )
ggsave(file.path(fig_dir, "color", "naacp_wartime_adjusted.pdf"), plot = p_color, width = 5, height = 3.5)

p_bw <- p_color +
  scale_fill_manual(values = c("No" = "#555", "Yes" = "#999"))
ggsave(file.path(fig_dir, "bw", "naacp_wartime_adjusted.pdf"), plot = p_bw, width = 5, height = 3.5)
